Reference 


NBS 

Publi- 

cations 


'AirWJii— 2-s-ferffETr — 


NAT'L  INST.  OF  STANO  & TECH  R.I.C. 


A11104  ESTfibE 


NBSIR  78-1509 


An  Algorithm  and  Basic 
Computer  Program  for 
Calculating  Simple  Coal 
Gasification  Equilibria 


...  / 

William  S.  Horton 


Chemical  Stability  and  Corrosion  Division 
Center  for  Materials  Science 
National  Measurement  Laboratory 
National  Bureau  of  Standards 
Washington,  D C.  20234 


NATIONAL  BUREAU  OF  STANDARDS 


a*mJ  Wof  Standards 

APR  17  197§ 


NBSIR  78-1509 

» *• 

AN  ALGORITHM  AND  BASIC 
COMPUTER  PROGRAM  FOR 
CALCULATING  SIMPLE  COAL 
GASIFICATION  EQUILIBRIA 


William  S.  Horton 


Chemical  Stability  and  Corrosion  Division 
Center  for  Materials  Science 
National  Measurement  Laboratory 
National  Bureau  of  Standards 
Washisngton,  D C.  20234 


August  1978 


MCl&UaU, 


OSc 


U.S.  DEPARTMENT  OF  COMMERCE,  Juanita  M.  Kreps,  Secretary 


Dr.  Sidney  Harman,  Under  Secretary 

Jordan  J.  Baruch,  Assistant  Secretary  for  Science  and  Technology 
NATIONAL  BUREAU  OF  STANDARDS,  Ernest  Ambler.  Director 


1 


I 

I 


' 


TABLE  OF  CONTENTS 


Page 

Introduction  1 

Thermodynamic  Basis  and  Assumptions  4 

Thermodynamic  and  Mathematical  Treatment  7 

General  Statement  of  the  Algorithm  17 

Application  to  a Simplified  Coal  Gasification  System  25 

Thermodynamic  Data  30 

Detailed  Algorithm  32 

Discussion  of  the  Listing  45 

Information  for  Using  the  Program  47 

The  Phase  Rule  50 

Examples  55 

Appendices  53 

1.  Symbols  - Text  53 

2.  Symbols  - Listing  55 

3.  Listing  57 


References 


75 


An  Algorithm  and  BASIC  Computer  Program  for  Calculating  Simple 
Coal  Gasification  Equilibria 

William  S.  Horton 
Center  for  Materials  Science 
National  Bureau  of  Standards 
Washington,  DC  20234 

INTRODUCTION 

Designing  plants  for  the  gasification  of  coal  is  being  aided  by 
testing  materials  expected  to  be  capable  of  withstanding  the  required 
conditions.  A variety  of  physical  and  chemical  tests  are  performed 
under  these  possible  conditions  and  in  the  same  gaseous  environments. 
Generally,  however,  to  fully  understand  the  observed  effects,  tests  are 
made  in  a variety  of  environments  related  to  those  of  the  gasification 
plants  and  at  a variety  of  conditions  of  temperature  and  pressure.  The 
gas  mixtures  used  as  the  environments  are  generally  not  at  chemical 
equilibrium,  and  depending  upon  the  conditions  and  the  possible  presence 
of  catalyzing  substances,  chemical  reactions  will  take  place  [2]  at  some 
unpredictable  speed.  The  tester  is  concerned  about  the  true  nature  of 
the  environment,  but  this  can  only  be  found  by  carefully  arranged  chemical 
analyses.  On  the  other  hand  it  is  of  value  to  know  the  maximum  change 
that  can  occur  during  the  test.  This  can  be  determined  if  the  equi- 
librium composition  for  the  mixture  can  be  calculated. 
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This  report  is  concerned  with  a simple  procedure  for  calculating 
the  equilibrium  composition  of  such  testing  mixtures.  The  aims  in 
devising  this  procedure  were  to  provide  an  interactive  computer  program 
that  was  easy  to  use,  not  inherently  complex,  that  exploited  the  relative 
simplicity  of  the  BASIC  programming  language,  and  was  confined  to  the 
coal  gasification  problem  alone.  Such  a program  would  be  usable  by  non- 
programmers, who  merely  need  to  sign  on  an  interactive  timesharing 
terminal  and  type  answers  to  a few  questions  about  the  testing  conditions. 
There  are  batch  programs  already  available  for  solving  simultaneous 
chemical  equilibria,  e.g.,  CEC71  [8],  however  these  have  many  disadvantages 
for  the  problem  at  hand.  They  are  often  complex,  handle  many  more 
problems  than  are  needed,  require  extensive  storage,  have  excessive 
print-out,  and  require  some  learning  prior  to  use.  There  is  an  interest- 
ing algorithm  [14]  with  a computer  program  written  in  ALGOL,  which 
treats  multiphase  chemical  equilibria  involved  in  coal  combustion.  This 
includes  formation  of  slag. 

The  development  of  a program  fulfilling  these  aims  is  presented 
here.  Included  are  the  mathematical  and  chemical  thermodynamic  treatment. 
Instructions  for  use  and  a listing  are  provided.  Also  included  are  some 
results  obtained  using  the  program  which  illustrate  the  use  of  the  phase 
rule  to  predict  the  potential  condensation  of  solid  carbon  during  a 
materials  test  under  coal  gasification  related  conditions. 

The  program,  called  COLGAS,  is  capable  of  calculating  equilibrium 
gas  compositions  containing  CH^,  CO,  CO2,  and  h^O  as  well  as  the 
quantity  of  solid  carbon  and/or  liquid  water  when  these  are  present. 
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Constant  temperature  is  assumed,  but  the  operator  may  choose  to  hold 
either  pressure  or  volume  constant. 

The  exposition  may  appear  more  detailed  than  is  warranted  for  the 
simplicity  of  the  problem  to  which  the  development  is  applied.  The 
author  found  this  to  be  helpful  in  seeing  the  relationship  of  the  parts 
and  the  exact  nature  of  the  assumptions  used.  The  detail  has  been 
retained  in  the  report  with  the  hope  that  some  readers  may  benefit 
thereby. 

Wherever  possible  the  recommendations  of  The  International  Union  of 
Pure  and  Applied  Chemistry  (IUPAC)  for  physicochemical  symbols  and  units  [3] 
have  been  followed.  This  includes  the  SI  units.  However,  the  atmo- 
sphere, liter,  and  degree  Celsius  have  been  made  usable.  The  SI  units 
used  are  the  megapascal  (MPa)  and  the  kelvin  (K).  Because  of  the  desire 
to  keep  symbolism  the  same  in  the  text  and  in  the  computer  program 
wherever  possible,  occasionally  the  IUPAC  symbol  has  not  been  used.  For 
example,  k is  used  for  equilibrium  constant  rather  than  K. 
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THERMODYNAMIC  BASIS  AND  ASSUMPTIONS 

The  literature  dealing  with  calculation  of  chemical  equilibrium 
compositions  is  abundant,  and  most  of  it  up  to  1970  is  given  by  van 
Zeggeren  and  Storey  [9].  Other  useful  reviews  are  by  Zeleznik  and 
Gordon  [13]  and  Klein  [18].  The  two  known  approaches  depend  upon  either 
solving  simultaneous  equilibrium  constant  equations  or  minimizing  the 
Gibbs  energy  of  the  system  (for  a constant  pressure  problem).  Both 
include  the  constraints  for  conserving  the  chemical  elements.  The  first 
of  these  methods  is  rather  old,  but  Brinkley  [10]  is  credited  with 
developing  a generalized  scheme  for  a systematic  treatment  using  equilibrium 
constants.  This  procedure  is  amenable  to  the  use  of  computers.  The 
first  reference  to  minimizing  the  Gibbs  energy  as  a route  to  solving 
simultaneous  chemical  equilibria  was  probably  by  White,  Johnson,  and 
Dantzig  [11].  In  principal,  however,  these  two  approaches  are  equivalent. 
This  becomes  obvious  if  minimization  is  effected  by  setting  the  partial 
derivatives  with  respect  to  the  constituents  equal  to  zero.  The  result 
is  a set  of  equations  involving  the  logarithms  of  the  species  concen- 
trations. If  logarithms  are  taken  of  the  Brinkley  equilibrium  constant 
equations  a set  of  similar  equations  is  obtained.  With  a certain 
choice  of  "constituents"  in  the  Brinkley  sense,  these  sets  of  equations 
are  identical.  A real  difference  in  the  methods  can  occur  if  mini- 
mization is  effected  by  other  methods  such  as  steepest  descent  or  linear 
programming  [11]. 

Minimization  of  the  Gibbs  energy  by  use  of  the  partial  derivatives 
has  at  least  one  advantage  over  the  Brinkley  method.  There  is  no  need 
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to  choose  substances  likely  to  be  present  in  larger  abundance  than 
others  and  to  treat  these  differently.  All  substances  in  the  gas  phase 
are  treated  alike  and  therefore  the  algorithm  can  be  a bit  less  complex. 
Furthermore,  the  addition  of  a new  substance  is  always  done  in  the  same 
way  and  does  not  lead  to  rewriting  the  equations  already  in  the  algorithm. 
Only  new  equations  are  added.  For  these  reasons  minimizing  the  Gibbs 
energy  was  the  chosen  method.  Minimization  is  accomplished  by  the  use 
of  Lagrangian  multipliers  and  setting  the  partial  derivatives  equal  to 
zero  followed  by  the  Newton-Raphson  method  of  solution.  Minimization  of 
the  Helmholtz  energy  for  constant  volume  problems  leads  to  equations  only 
slightly  different  from  those  for  constant  pressure  and  to  minor  additions 
to  the  program. 

For  the  testing  of  coal  gasification  plant  materials  a reference 
environment  has  been  chosen  [12]  and  is  shown  in  Table  1. 


Table  1.  Standard  Gas  Environment 
Volume  Percent 


CH4 

CO 

co2 

H2 

H20 

h2s 

NH3 

5 

18 

12 

24 

39.5 

0.5 

1 

The  presence  of  H^S  can  lead  to  S02  and  the  NH^  can  give  N2<  How- 
ever, all  four  of  these  will  be  in  relatively  low  concentration  and  it 
was  decided  to  consider  the  first  five  species  only.  If  desired,  the 
concentrations  of  the  others  can  be  simply  calculated  once  these  five 
are  known,  as  well  as  the  amount  of  sulfur  and/or  nitrogen. 


6 


All  gases  are  assumed  to  be  ideal.  Condensed  water  and  carbon  are 
assumed  to  be  pure.  For  example,  the  solubility  of  carbon  dioxide  in 
water  is  ignored.  For  constant  volume  problems,  the  volume  of  condensed 
species  is  ignored  in  comparison  to  the  gas  volumes.  That  this  is 
generally  acceptable  is  evident  from  the  fact  that  when  90%  of  the 
original  gas  has  been  equilibrated  as  condensed  species,  these  latter 
contribute  only  about  0.5%  to  the  total  volume  in  the  worst  case. 
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THERMODYNAMIC  AND  MATHEMATICAL  TREATMENT 


Assume  a system  of  s substances,  g of  which  are  gases,  at  pressure 
P,  with  volume  V,  and  thermodynamic  temperature  T.  If  the  system  is  not 
at  chemical  equilibrium,  chemical  reactions  will  adjust  the  number  of 
moles  of  each  substance,  n-,  until  equilibrium  is  attained.  At  this 
point  the  Gibbs  energy,  G,  will  have  become  a minimum  if  the  pressure 
and  temperature  are  kept  constant  (ref  1,  p 209).  If,  instead,  the 
volume  and  temperature  are  kept  constant,  the  Helmholtz  energy,  A,  will 
have  become  a minimum,  (ref  1,  p 208).  The  s-g  > 0 substances  are 
condensed  phases  and  are  assumed  to  be  pure  and  immiscible  with  one 
another.  The  gases  are  assumed  to  be  ideal  and  insoluble  in  the  condensed 
phases. 

The  Gibbs  energy  may  be  expressed  as  the  sum  of  the  contributions 
from  each  substance. 


The  p.  are  Gibbs  energies  per  mole  and  are  known  as  the  partial  molar 
Gibbs  energies  or  the  chemical  potentials  (ref  1,  pp  214,215).  For  the 
g gases  in  an  ideal  mixture  the  chemical  potentials  are  a function  of 
the  partial  pressures,  p • (ref  1,  p 261). 


s 


(1) 


p.  = p?  + RT  In  p. 


(2) 
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The  p°  are  constants  at  fixed  temperature  and  represent  the  Gibbs 
energy  of  one  mole  of  pure  ideal  gas  at  0.101325  megapascals  (one  atmo- 
sphere) pressure.  For  pure  condensed  phases  the  p.  are  merely  the  p° 
for  these  substances.  Consequently,  equation  (1)  can  now  be  rewritten. 

s g 

G = H n.p°  + RT  I n.ln  p.  (3) 

1 1 1 1 1 1 

Note  that  the  first  g subscripts  have  been  assigned  to  the  gases.  The  partial 
pressures  are  equal  to  the  product  of  the  mole  fractions,  n^./n^,  anc*  the 
total  pressure. 

p.  = (n./n  ) P 

pi  i g' 

g 

n = \ ' n. 

g L*t  i 

1 

Thus,  equation  (3)  may  be  rewritten. 

s g 

G = X"iM°  + RT  XX  ln  (niP/V 
1 1 

The  p.  may  be  evaluated  by  considering  compounds.  For  example,  the 
compound  AxB^  may  be  made,  in  principle,  from  its  elements,  A and  B, 
according  to  the  following  chemical  reaction. 

xA  + yB  ->  A B 
J x y 


(4) 


(5) 


I 
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The  Gibbs  energy  change  for  the  reversible  formation  of  one  mole  of 
A^By  is  given  by  AG°  with  all  the  substances  in  the  standard  states  (ref 
1,  p 283). 

AGf  = mA  B ’ x mA  " y mB 
x y 

But  for  an  equilibrium  among  these  three  substances,  the  AG^  is  related 
to  the  equilibrium  constant  (ref  1,  p 283). 

AG°  = - RT  In  k. 


He 


re  k.  is  the  equilibrium  constant  for  formation  of  one  mole  of  sub- 


stance and  R is  the  molar  gas  constant.  Consequently,  one  can 

express  the  chemical  potential  of  a compound  in  terms  of  the  absolute 
chemical  potentials  of  its  elements  and  the  equilibrium  constant  of 
formati on. 

mA  B = X mA  + y PB  " RT  1n  kf 
x y 

For  the  general  case  let  k.  represent  the  equilibrium  constant  of 
formation  for  the  i - th  substance. 


'i  S Vi  J 
j=l 


RT  In  k. 


(7) 


The  summation  is  over  all  e elements  present  in  the  system,  with 

v.  . = 0 for  any  element,  j,  not  present  in  substance  i.  Furthermore, 
^ J 

this  expression  may  be  used  for  substances  which  are  themselves  chem- 
ical elements  by  recognizing  that  for  these  In  k.  = 0.  Using  these 
ideas,  the  first  term  in  equation  (6)  may  be  rewritten. 


RT  In  k. 

i 


i=l 


j=l 
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It  can  be  shown  that  the  right  hand  side  here  can  be  arranged  as 
fol lows. 

h Vi = i;nivij-RT  Eniinki  («) 

1 j=l  i=l  i=l 

Because  the  elements  are  neither  destroyed  nor  created  by  shifts  among 
the  n.j,  S n_j  v — - E j , the  number  of  gram  atoms  of  element  j originally 
introduced  into  the  system.  Consequently,  the  first  term  on  the  right- 
hand  side  of  equation  (8)  is  given  by 

H E.  = G° 
j=l  J J 

where  G°  is  the  absolute  Gibbs  energy  due  to  the  elements,  a constant. 
Therefore,  the  sum  of  Gibbs  energies  for  all  the  substances  in  standard 
states  is  given  by  the  following. 

s s 

£ n.q?  = G°  - RT  En.  In  k.  (9) 

1 1 1 1 1 1 

This  may  be  used  in  (6). 

s g 

G = G°  - RT  S n.  In  ^ + RT  S n.  In  (n.  P/ng)  (10) 

The  right-hand  member  may  be  rearranged  to  give  a new  function, 
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9 s 

rG  = TT  = £ "i  ln  {ni  P/V  ' ? ni  ln  ki  <") 

I 3 1 

Now,  because  G°,  R,  and  T are  constants,  minimizing  will  minimize  G. 

Turning  at  this  point  to  the  problem  of  equilibrium  at  constant 
volume  and  temperature,  use  is  made  of  the  relationship  between  G and  A 
(ref  1,  p 202). 

A = G - PV  (12) 

Use  of  equation  (10)  gives  an  expression  which  will  lead  to  one  that  is 
si  mi  1 ar  to  that  in  (11). 

s g 

A = G°  - RT  y n.  In  k.  + RT  y n.  In  (n.  P/n  ) - PV  (13) 

i i Ar*  i i q 

1 1 y 

Assuming  that  the  volume  of  the  condensed  phases  is  negligible  compared 

to  that  of  the  gases  permits  using  the  (ideal)  gas  volume  for  V. 

PV  = ng  RT  (14) 

At  the  same  time  it  is  possible  to  replace  the  variable  P. 

s g 

A = G°  - RT  E n.  ln  k.  + RT  Vn.  In  (n.RT/V)  - n RT  (15) 

l i ' ii  g 

1 1 

In  view  of  equation  (5)  the  last  term  may  be  put  inside  the  gas  sum- 
mation and  the  expression  rearranged  to  define  the  function  = 

9 S 

rA  = ^1  = E n{  [In  (n.  RT/VJ-1]-  L n.  ln  ki  (16) 
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Minimization  of  V ^ also  minimizes  A. 

The  problem  of  equilibration  has  now  been  shown  to  require  min- 
imization of  rG  or  PA  respectively  for  constant  pressure  or  constant 
volume.  However,  these  minimizations  with  respect  to  the  n.  must  be 
carried  out  subject  to  constraints:  1.  The  chemical  elements  are 
conserved  as  previously  noted,  and,  2.  No  n-  may  be  less  than  zero. 

The  constraints  may  be  expressed  as  follows. 

s 

L = E n,  ,,  - E,  = 0 j = 1(1 )e  (17) 

J 1 1 1J  J 

n.  > 0 i = 1 (1 )s  (18) 


There  are  e functions,  L . , each  of  which  expresses  the  conservation  of 

element  j such  that  the  E-  equals  the  number  of  gram  atoms  of  that 

1 

element  in  the  system.  The  E.  are  determined  from  the  initial  com- 
position. As  before,  . . is  the  subscript  for  the  element  j in  sub- 

J 

stance  i.  Constraint  (18)  is  handled  in  this  development  by  replac- 
ing negative  estimates  of  gaseous  constituents  by  small  positive  num- 
bers. A negative  n.  for  a condensed  phase  assumed  to  be  present  indi- 
cates the  absence  of  that  phase  and  is  retained  as  an  indicator  of  that 
fact. 


The  method  of  minimization  adopted  here  is  the  classic  one  using 

Lagrangian  multipliers,  A.,  [6].  To  minimize  V form  the  function  F by 

J 

multiplying  each  of  the  L.  by  A.  and  adding  to  V. 

J J 
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Fr  = Z n.  In  (n.  P/n  ) - Z n.  In  k.  + Z A- 
b 1 1 9 -|  1 1 i J 


1 


1 


In, 


F.  = Z n.  In  (n.  RT/V)  - Z n.  In  k.  + Z A. 

" 1 1 l -jiT-iJ 


1 


, 'I  vifEj 


£ n.  v.  .-E. 
1 T 1J  J 


(19) 


(20) 


J 


Minimization  of  F^  will  be  treated  first  in  detail.  It  will  be  seen 
later  that  a few  simple  changes  in  the  equations  will  suffice  to  minimize 

FA- 

Setting  the  appropriate  partial  derivatives  of  F^  equal  to  zero 
results  in  three  sets  of  equations,  one  set  for  the  gases,  one  for  the 
condensed  phases,  and  one  for  the  elements. 

e 


f . = In  n.  - In  n.  + 


g ln  (PAi)  + £ «ij  = 0 
j=l 


i=l(l)g  (2 1) 


f . = In  k.  + 


A • vi  • - — 0 


. + V A . „ 

1 2-  J 1 J 

j=l 


i=(g+i)(i)s  (22) 


f_-  = V n.  v - -•  " E, 
i=l 


t Z i vij  j 


j=l (1 )e  (23) 


It  can  be  seen  that  there  are  s + e equations,  sufficient  to  completely 

determine  the  s n.  and  the  e A..  The  variable  n is  a defined  function 
i J 9 

which  is  calculated  when  the  n.  have  been  estimated. 

Equations  (21)  are  nonlinear  and  an  iterative  procedure,  Newton- 
Raphson  [7],  is  used.  In  summary,  this  method  iterates  from  initial 
estimates  of  the  n-  and  the  A.-  to  a final  set  which  satisfies  some 

J 

stopping  criterion.  The  iterations  may  be  expressed  by  use  of  a matrix 
equation  [13]. 
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(24) 


Here  x and  f are  column  matrices  (vectors)  representing  the  unknowns  and 
the  functions  (21)  through  (23).  Subscript  t indicates  that  the  values 
used  are  from  the  t-th  trial.  J is  the  Jacobian  of  the  f's,  a square 
matrix  with  elements  equal  to  the  partial  derivatives  of  the  functions, 
f,  with  respect  to  the  and  the  A . . 

There  are  a variety  of  types  of  elements  in  the  Jacobian.  Equations 
(21)  yield  the  following  entries. 


3f . 

l 


n 


i 


i=l(l)g  (25) 


3f 


3n 


1 = - 1/n 


|i=KDg 

jk=l(l)s  (26) 

' k^i 


V 

i j 


Equations  (22)  yield  the  following. 


ji=KDg 

/j=l(l)e  (27) 


3f . 

l 

3n, 

k 


= 0 


(i=g+l(l)s 
| k=l (1 )s 


(28) 


3f . 

3A  • 
J 


v 

i j 


j i=g+l(l)s 
j j=l(l)e 


(29) 


Equations  (23)  yield  the  following. 
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8n . 


i 


jj=l(l)e 
Vij  f i=l(l  )s 


(30) 


8f  - 


k 


j— 1 ( 1 )e 
k=l (1 )e 


(31) 


Expression  (25)  gives  nonzero  values  to  the  first  g Jacobian  elements 
along  the  diagonal.  Expressions  (27)  and  (29)  may  be  combined. 


Expressions  (26),  (30),  and  (32)  show  that  these  elements  are  symmet- 
rically disposed  with  respect  to  the  diagonal.  Since  all  other  elements 
are  zero,  J is  symmetric.  This  is  an  advantage  for  some  methods  of 
calculating  the  inverse. 

These  results  are  easily  modified  for  equilibria  at  constant  volume. 
The  functions  (21)  are  replaced  by  (33),  but  (22)  and  (23)  are  retained. 


i=l (1 )s 
j=l(l)e 


(32) 


e 


(33) 


Of  the  Jacobian  elements,  those  in  (25)  and  (26)  are  changed,  replaced 
by  (34)  and  (35). 
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3f  • 

57T  = 1/n- 
9n^.  i 

i=KDg  (34) 

af . 

--  = 0 
3nk 

i=Ki)g 

k=l (1 )s  (35) 

k*i 
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GENERAL  STATEMENT  OF  THE  ALGORITHM 


In  principle  there  is  no  difficulty  obtaining  equilibrium  composi- 
tions from  the  development  of  the  previous  section.  A problem  is  defined 
by  giving  values  for  the  constant  temperature  and  the  constant  pressure 
(or  volume)  and  for  the  composition  of  the  test  gas  in  moles  of  each 


chemical  element  may  be  calculated.  From  the  formulae  of  the  substances 


be  obtained  from  the  literature.  One  needs  only  initial  estimates  of 
the  n . and  the  A. . in  order  to  proceed. 

1.  Choose  initial  estimates. 

2.  Evaluate  the  functions,  equations  (21)  to  (23),  using  the 
initial  estimates. 

3.  Calculate  the  elements  of  the  Jacobian,  equations  (25)  to  (32). 

4.  Invert  the  Jacobian. 


5.  Calculate  the  new  estimates  from  equation  (24). 

6.  Check  the  estimates  against  the  stopping  criterion. 

7.  Depending  upon  the  result  of  step  6,  either  stop  or  repeat 
steps  2 through  6 using  the  new  estimates  for  evaluation 
purposes. 

* 

In  the  later  section  about  using  the  program  the  instructions  ask  for 

mole  fractions  or  percentages  although  these  data  are  treated  as  moles. 

Experience  with  COLGAS  has  shown  that  when  n is  much  less  or  much 

9 

greater  than  one,  difficulty  with  over-  or  underflow  may  occur. 


substance. 


the  v..  are  also  determined. 


The  equilibrium  constants  of  formation  may 
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As  straightforward  as  this  procedure  is,  there  are  a number  of 
practical  points  associated  with  these  steps  which  are  worth  commenting 
on.  However,  the  handling  of  the  practical  problems  depends  upon  the 
general  aims  of  the  program.  Therefore,  although  some  general  remarks 
will  be  made,  it  must  be  understood  that  the  particular  procedure 
adopted  in  this  report  is  based  upon  the  restricted  nature  of  the 
problems  being  solved. 

Choosing  the  Estimates.  The  program  to  be  developed  was  to  give  a "one 
shot"  solution.  That  is,  a single  run  would  produce  only  a single 
answer.  A different  temperature  and/or  a different  pressure  (or  volume) 
or  test  gas  composition  will  require  a new  run.  Consequently,  the 
program  has  no  prior  experience  to  build  upon  as  is  used  in  some  programs 
[8,  14].  It  was  decided  to  use  the  test  gas  composition  as  the  initial 
estimates  for  the  substances,  and  this  has  turned  out  to  be  satisfactory. 
For  each  A.,  zero  has  been  a satisfactory  initial  estimate. 

There  are  two  versions  of  the  program  in  the  same  listing.  The 
commonly  used  version  has  a short  printout.  Either  the  short  or  the 
long  printout  is  obtained  by  an  appropriate  response  to  the  initial 
query  at  the  terminal.  The  long  version  permits  the  operator  to  suggest 
which  of  the  four  combinations  of  condensed  species  situations  is  believed 
to  be  correct.  It  then  allows  input  of  a set  of  initial  estimates 
different  from  that  used  automatically  by  the  program.  These  alternatives 
have  been  provided  so  that  on  occasions  of  non-convergence  an  operator 
may  try  to  redo  the  iterations  with  a better  starting  situation. 
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Evaluating  the  Functions.  Equations  (21)  to  (23)  present  no  problems 
beyond  obtaining  values  for  the  equilibrium  constants  of  formation  for 
compounds  in  the  system.  (The  total  number  of  moles  of  gases,  n , is 
evaluated  anew  at  each  iteration  using  equation  (5)).  The  equilibrium 
constants  used  in  the  program,  except  for  liquid  water,  come  from  the 
JANAF  Tables  [15]  and  are  discussed  in  a separate  section. 

Calculating  the  Elements  of  the  Jacobian.  Examination  of  equations  (25) 
to  (32)  shows  that  the  only  elements  varying  in  value  from  iteration  to 
iteration  are  the  elements  of  the  upper  left  5x5  submatrix.  These  will 
require  recalculation  during  the  iteration.  All  other  elements  are 
constant  integers  and  need  not  be  re-evaluated.  However,  if  the  presence 
or  absence  of  a condensed  phase  turns  out  to  be  different  than  assumed, 
the  presence  or  absence  of  the  appropriate  constant  elements  must  be 
adjusted. 

Inverting  the  Jacobian.  There  is  an  extensive  literature  on  the  inversion 
of  matrices,  particularly  symmetric  ones  [16].  It  is  tempting  to  seek 
a direct  analytical  solution  with  perhaps  the  possibility  of  avoiding 
some  round-off  error  in  the  computer.  In  fact,  at  constant  volume  and 
for  the  case  where  solid  carbon  is  present,  but  liquid  water  is  not,  it 
is  relatively  simple  to  do  this  using  either  an  LP  or  an  LDR  decomposition 
of  the  Jacobian  [17].  There  is  a minor  advantage  to  this  beyond  the 
possibility  of  reducing  the  round-off  error.  Some  of  the  expressions 
for  the  n . contain  functions  of  them  which  may  be  replaced  by  the  Ej 
constants  for  hydrogen  and  oxygen.  This  should  improve  the  accuracy  of 
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the  estimates  and  perhaps  reduce  the  number  of  required  iterations. 
However,  when  the  condensed  phase  situation  is  reversed,  some  of  the 
expressions  for  the  estimates  become  quite  complicated  and  hardly  worth 
evaluating  directly.  Although  not  essential,  it  is  desirable  to  handle 
all  varieties  of  condensed  phase  situations  in  the  same  way  for  sim- 
plicity of  programming.  This  is  particularly  true  when  using  the  BASIC 
programming  language  because  of  the  one-statement  matrix  inversion 
procedure  provided.  A direct  analytical  solution  was,  therefore,  not 
used. 

Calculating  the  New  Estimates.  As  is  well  known,  the  Newton-Raphson 
procedure  often  overcorrects  in  the  early  stages,  particularly  if  the 
initial  estimates  are  somewhat  far  from  the  mark.  This  results  in 
divergence.  To  avoid  this  a common  procedure  is  to  modify  the  correct- 
ion term  in  the  estimating  equation,  (24),  by  a weighting  factor,  k, 
less  than  unity. 


xt+i  = xt 


kt  V ft 


(36) 


In  many  nonlinear  equation  problems  it  is  adequate  to  use  k = 0.5  in  the 
early  stages  and  to  change  it  to  unity  later  on.  In  running  this 
program  it  was  found  that  this  was  not  always  sufficient.  This  diffi- 
culty was  overcome  by  making  k a variable  function  of  t,  the  number  of 
the  trial.  A sinusoidal  function  of  t sets  k equal  to  about  0.01  for  t 
= 1,  to  0.5  for  t = 8,  and  to  1.0  for  t = 16.  The  value  of  k remains  at 
unity  for  any  further  iterations.  This  is  done  only  for  the  gaseous 
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constituents.  Experience  with  the  algorithm  has  shown  that  the  estimates 
of  the  Langrangian  multipliers  converge  easily  with  k = 1.0  at  all 
stages. 

Sometimes  new  estimates  of  the  gaseous  moles  become  negative.  This 
not  only  violates  the  constraint,  equation  (18),  but  it  cannot  be  tol- 
erated when  evaluating  the  functions,  f . , which  contain  In  n..  Trans- 

2 

formations  may  avoid  this  problem,  y = n [4]  or  y = In  n [5],  but  intro- 
duce numerous  back  and  forth  transformation  steps  in  the  algorithm  which 
were  thought  to  be  undesirable.  The  algorithm  reported  here  replaces 
such  numbers  with  a very  small  number,  originally  1 x 10  ^ . It  was 
found  that  this  occasionally  led  to  difficulty  when  the  size  of  the 
correction  was  larger  than  this  figure.  A loop  of  repeated  values  could 
result  with  continually  negative  estimates  for  the  same  n-  leading  to  a 
lack  of  convergence.  This  difficulty  was  partially  overcome  by  replacing 
the  7 in  the  exponent  by  a random  number  distributed  uniformly  between  7 
and  10  inclusive.  This  choice  occasionally  leads  to  a new  difficulty 
with  some  processors.  During  the  matrix  inversion  step,  usually  a 
Gauss-Jordan  routine,  the  product  of  the  five  (1/n.j-l/n  ) is  evaluated. 

If  for  a given  iteration  a number  of  n^  estimates  become  negative  the  n^ 
product  may  cause  an  overflow.  This  has  happened  with  a compiler  for 
which  the  exponent  of  numerical  constants  is  restrained  to  be  within  ±38. 
For  a processor  with  a larger  range,  say  greater  than  ±55,  this  difficulty 
may  not  arise. 
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Checking  the  Estimates  for  Stopping.  The  criterion  used  here  for 

stopping  depends  upon  the  magnitude  of  the  unweighted  correction  to  the 

n.  and  the  A..  The  vector  c.  is  the  correction  and  it  comes  from  (24). 
i J t 

ct  = a1  ft 

The  iterations  are  considered  complete  when 

ct  -5 

— — < 5 x 10  3 (38) 

nt 

for  all  n.j  and  Aj.  This  particular  stopping  criterion  has  been  crit- 
icized, (ref  18,  p 520)  as  not  leading  necessarily  to  the  same  solution 
as  a criterion  based  upon  the  derivatives:  f^  being  sufficiently  close 
to  zero.  However,  because  the  interest  here  is  in  the  gas  composition 
the  criterion  given  above  is  preferable.  To  assure  adequate  minimi- 
zation, the  pertinent  independent  chemical  equilibria  constants  are 
checked  for  agreement  when  calculated  from  Gibbs  energies  of  formation 
and  from  the  calculated  mole  fractions.  According  to  Brinkley  [23] 
there  are  at  most  R = N - C reactions  for  N chemical  species  and  C com- 
ponents. In  this  discussion  C = 3 and  N equals  5,  6,  or  7 depending 
upon  the  case.  For  case  IV,  all  gas,  (see  below),  N = 5 and  R = 2. 

The  reactions  chosen  were  the  following. 

A:  C02  + H?  = CO  + H20 

B:  CO  + 3H2  = CH4  + H20 

When  C(s)  is  present,  reaction  C is  included,  and  for  liquid  water 


reaction  D. 
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C:  CG2  + C(s)  = 2C0 

D:  H20(1)  = H20(g) 

Terminating  the  Iterations.  When  the  stopping  criterion,  equation  (38) 

has  been  satisfied  for  all  n.  and  A.  one  set  of  iterations  is  complete. 

J 

The  problem  may  not  be  solved,  however.  In  order  to  perform  the  calcu- 
lations a choice  had  to  be  made  among  the  possible  combinations  of 
condensed  species  situations  as  shown  in  Table  2. 

Table  2.  Condensed  Phase  Cases 
H20(1) 

YES  NO 

YES  1-8  II  - 7 

C(s) 

NO  III  - 9 IV  - 8 

Here  the  yes-no  tabulations  refer  to  the  presence  or  absence  of  the 
phase  at  the  top  or  at  the  left.  The  Roman  numerals  are  used  to  identify 
the  cases,  and  the  Arabic  numerals  show  the  number  of  functions,  f., 
and  correspond  to  the  order  of  the  relevant  Jacobian. 

Solving  the  seven-equation  system  correspond!' ng  to  case  II,  say, 
has  assumed  the  presence  of  solid  carbon  and  the  absence  of  liquid 
water.  However,  this  assumption  may  not  be  correct.  Therefore,  after 
going  through  steps  1 to  7 and  deciding  to  stop,  it  is  necessary  to 
check  the  assumption  of  case  II.  If  liquid  water  is  really  present 
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and/or  solid  carbon  is  really  absent  as  suggested  by  the  final  n. , then 
the  appropriate  case  must  be  chosen  and  the  iterations  redone.  There 
are,  therefore,  two  hierarchies  of  iterations,  one  on  the  n.  and  one  on 

J 

the  cases.  When  the  final  n.  are  found  to  be  consistent  with  the  case 
assumed,  the  calculation  is  finally  finished.  Details  of  checking  case 
consistency  are  covered  in  the  section  on  DETAILS  OF  THE  ALGORITHM. 
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APPLICATION  TO  A SIMPLIFIED  COAL  GASIFICATION  SYSTEM 

The  work  shown  in  this  report  was  undertaken  to  provide  results 
with  atmospheres  used  for  testing  materials  to  be  used  in  coal  gasification 
plants.  For  simplicity  the  problem  was  limited  to  one  in  which  an 
initial  gas  mixture  was  made  from  two  or  more  of  the  following:  CH^,  CO, 
CC^,  and  F^O  at  some  temperature,  T.  The  test  may  be  carried  out 
dynamically  in  which  case  constant  pressure  is  assumed.  A static  test 
implies  constant  volume.  Constant  pressure  will  be  assumed  here.  From 
what  has  gone  before,  extension  to  constant  volume  will  be  obvious. 

The  first  step  is  to  set  down  the  pertinent  functions,  according  to 
equations  (21)  through  (23)  of  the  second  section,  but  it  is  necessary 
to  decide  upon  allowable  condensed  states.  Because  conditions  of  test 
will  not  permit  the  existence  of  ice,  only  liquid  water  need  be  considered. 
The  conditions  will  be  above  the  critical  temperatures  of  the  other  four 
substances  so  no  other  solid  or  liquid  phase  for  these  need  be  considered. 
Flowever,  decomposition  of  CH^  and  disproportionation  of  CO  can  produce 
solid  carbon.  Consequently,  there  are  two  condensed  phases  to  consider, 
liquid  water  and  solid  carbon.  Examination  of  the  simultaneous  equations 
in  the  light  of  this  shows  that  n , the  number  of  moles  of  solid  carbon, 
appears  only  once,  and  that  in  the  equations  (23)  corresponding  to  that 
element.  Allowing  indices  one  through  five  to  refer  to  the  gases  in  the 
order  mentioned  above 


nn  + n0  + n0  + n_  - 
12  3c 


C = 0 


(39) 
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Here  C is  the  E.  for  the  element  carbon.  It  is  unnecessary,  therefore, 
J 

to  solve  for  nc  iteratively.  Instead  the  equation  for  the  constraint 
due  to  the  element  carbon  will  be  deleted.  When  the  other  n.  have  all 
been  determined,  nc  may  be  calculated  from  (39).  For  case  II  (solid 
carbon  only)  the  functions  to  be  used  are  the  following. 


f , = In  n,  - In  n - In  kn  + 4A, 

1 1 g 1 h 


f0  - In  n0  - In  n - In  k„  + A 
2 2 g 2 o 


f,  - In  n.  - In  n - In  k0  + 2A 
3 3 g 3 o 


f4  = 1n  n4  " ln  n g ~ 1n  k4  + 2*h 


(40) 


fr  = In  nc  ~ In  n * In  L + 2A,  + A 
5 5 g 5 h o 


f6  = 4n]  + 2n4  + 2n5  “ h 


f7  ~ n2  + 2n3  + n5  ' 0 


Here  h and  o are  the  E.  for  hydrogen  and  oxygen. 

Equations  (25)  to  (31)  of  the  second  section  lead  simply  to  the 


Jacobian. 
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J = 


1 _ J_ 

n,  n 
1 9 


1 1 


n2  ng 


1 1 


n0  n 
3 g 


1 1 


n-  n 

4 g 


1 1 


nc  n 

5 g 


-1 


0 |(41) 


The  functions  of  equation  (40)  and  the  inverse  of  the  Jacobian  (41) 
are  the  needed  entries  for  the  right  hand  side  of  equation  (24)  of  the 
second  section. 
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The  development  of  case  II  gives  the  smallest  set  of  equations  to 
be  solved  by  iteration  and  give  rise  to  the  "basic"  Jacobian.  This  is  a 
7x7  upper  left  submatrix  common  to  the  cases.  Consequently,  the  elements 
in  it  may  be  defined  without  regard  to  case. 

The  development  here  was  for  a particular  case.  The  one  with 
condensation,  Case  I,  with  both  condensed  phases  present,  follows  by 
adding  terms  for  liquid  water  to  fg  and  f^:  2n^  and  n,  respectively.  A 
function  for  liquid  water  is  needed. 


The  pertinent  Jacobian  is  equation  (41)  with  another  row  and  column 
added.  The  row  is  given  by  jg. 


The  column  is  the  transpose  of  this. 

Case  III,  with  liquid  water  but  no  solid  carbon,  requires  fg 

and  n but  has  no  n . Because  there  is  no  n equation  (39)  becomes 
w c c 


This  function  must  be  included  as  well  as  a A . A new  row  and  column 

c 

are  added  to  the  Jacobian  consisting  of  ones  in  the  first  three  elements 
and  zeros  elsewhere.  This  is  the  largest  matrix  and  has  nine  rows  and 


(42) 


j8  - {0,0,0,0,0,2,1,01 


(43) 


(44) 


columns . 
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Case  IV,  with  neither  liquid  water  nor  solid  carbon,  is  derivable 
from  Case  III  by  elimination  of  the  eighth  row  and  column.  These  corr- 
espond to  liquid  water  and  A^.  Also  the  function,  fg,  is  deleted.  The 
problem  of  determining  equilibrium  compositions  in  the  system  is  then 
resolved  to  choosing  the  appropriate  form  of  the  Jacobian,  (41),  depend- 
ing upon  the  case.  Deciding  the  appropriate  case  for  any  given  problem 
without  any  prior  knowledge  is  difficult.  As  was  seen  in  the  previous 

y 

section,  the  algorithm  given  here  assumes  a starting  case  and  checks  the 
results  for  consistency  with  this  assumption. 
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THERMODYNAMIC  DATA 

For  the  five  gases,  CH^,  CO,  CO^,  H2,  H2O,  the  tabular  data  of  the 
JANAF  Tables  [15]  between  300  and  1400  K were  fitted  by  least  squares  to 
equations  of  the  form 

m 

In  k.  = a^  + a2  In  T + X]  ai+2  (45) 

Here,  as  stated  earlier,  k.  is  the  equilibrium  constant  of  formation  of 
the  i-th  compound  in  its  standard  state  and  from  its  elements  in  their 
standard  states.  The  constant  for  hydrogen  is,  of  course,  unity  so  that 
the  corresponding  In  k.  equals  zero.  The  value  of  m was  determined  by 
trials  with  m = 2(1)5.  The  smallest  m > 2 so  that  the  standard  error  of 
the  residuals  < .001  was  chosen.  The  form  of  this  equation  is  similar 
to  that  of  Barron,  Porter,  and  Hammond  [19]  but  the  latter  authors  used 
a fixed  highest  power  equal  to  five.  The  least  squares  coefficients 
were  determined  with  the  computer  language  OMNITAB  [20]  which  has  an 
accurate  least  squares  routine  [21]. 

Liquid  water  is  not  included  in  the  JANAF  Tables,  so  resort  was  had 
to  combining  the  JANAF  data  for  gaseous  H2O,  k^,  with  the  vapor  pressure, 

?s- 


1n  k£  = 1n  k5  " ln  ps 


(46) 


For  the  vapor  pressure  data,  Schmidt's  tabular  values  [24]  were  used  and 
fitted  by  least  squares  to  an  equation  of  the  following  form. 
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m 

In  ps  = b1  + b2/T  + E bi+2  T1  (47) 

i=l 

In  this  case  the  data  were  sufficiently  detailed  that  a tighter  criterion 
could  be  used  for  m.  The  smallest  m was  chosen  for  which  the  maximum 
residual  < 0.0005.  The  results,  In  for  the  gases  and  In  p&  for 
H20(£)  are  shown  in  the  table. 


Table  3.  Empirical  Coefficients  for  the  In  and  Ln  p$ 


In  k . 

ai 

a2 

a3xl03 

a^xlO^ 

9 

a^xlO3 

in13 

a^xlO 

:H4 

8212.35 

-0.816063 

-10.0801 

7.73804 

-3.06981 

5.08552 

:o 

13632.3 

1.8055 

- 2.3936 

0.372836 

— 

— 

,°2 

47326.7 

0.0777872 

- 0.280706 

0.0348747 

— 

— 

H2°(g) 

28786.6 

-0.697311 

- 1.45957 

0.83769 

-0. 180286 

— 

bl 

b2 

b3 

b^xlO4 

brXlO7 

5 

, nl  0 

b^xlO 

!"  ps 

55.5049 

-9430.74 

-0.159684 

3.07499 

-3.0549 

1.24325 

1 

The  maximum 

residual  for  any  of  the  nonzero  ln  k. 

to  be  used  is 

0.00186  for  ln  p 

. Program  trials  with  In 

k differing 

by  0.01  showed 

absolute  values 

of  changes 

in  the  computed 

n-j  < 0.0001. 

The  accuracy 

of  the  coefficients  in  the 

table  is  quite  adequate  for 

the  purpose. 
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DETAILED  ALGORITHM 

Master  Flowchart.  The  details  of  the  algorithm  will  be  described  in 
terms  of  flowcharts  for  some  sections  and  with  an  overall  chart  showing 
the  relationship  of  the  pieces.  It  will  be  noticed  that  although  some 
sections  of  the  program  could  be  subroutines,  they  are  not.  There  are 
no  "gosub"  and  "return"  statements.  There  is  no  advantage  to  having 
formal  subroutines  here.  Replacement  of  a section  is  no  more  difficult 
than  replacing  a subroutine  because  BASIC  statements  must  all  have  line 
numbers  which  must  be  adjusted  in  either  case.  Furthermore,  because  one 
cannot  fruitfully  use  dummy  variables,  one  of  the  advantages  of  a 
FORTRAN  type  subroutine  is  lost. 

The  master  flowchart  is  shown  as  Figure  1.  Most  of  the  items  are 
self-explanatory,  however,  some  remarks  about  others  will  be  useful. 

The  first  major  branch  point  reflects  the  impossibility  of  H2OO) 
existing  above  its  critical  temperature,  647.3  K,  and  below  that  at 
total  pressures  less  than  its  vapor  pressure.  This  is  a useful  distinction 
because  it  rules  out  two  of  the  four  cases,  I and  III.  The  choice  of  II 
and  IV  as  starting  cases  for  the  respective  branches  resulted  from  the 
experience  that  these  led  more  often  to  convergence  of  the  first  set  of 
iterations.  A longer  experience  may  suggest  alternate  choices.  Note 
that  the  variable  c-j  is  assigned  the  value  of  the  case  number  assumed  at 
the  start  of  a set  of  iterations.  This  is  required  for  the  case-consistency 
checking  section. 

The  Newton-Raphson  iterations  follow  in  either  case.  These  will  be 
treated  later  by  a separate  flowchart.  Upon  successful  convergence  of 
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the  iterations,  the  next  step  is  to  check  the  case  suggested  by  the 
results.  Upon  unsuccessful  convergence,  the  " lack-of-convergence 
output"  results.  The  latter  not  only  states  lack  of  convergence  but 
also  gives  the  case  number  that  had  been  assumed  and  the  last  set  of 
iterated  results.  These  outputs  are  useful  if  it  is  desired  to  try 
again  with  a different  set  of  initial  estimates. 

The  case-checking  section  will  be  presented  in  a detailed  flowchart. 
Essentially  it  computes  the  case  number,  c^,  suggested  by  the  success- 
fully converged  results  and  compares  this  with  c-j . If  they  are  equal 
then  a "successful  output"  results  which  gives  the  equilibrium  gas 
composition,  quantities  of  liquid  water  and/or  solid  carbon  if  either  or 
both  are  present,  and  the  changed  volume  (or  pressure).  If  c ^ does  not 
equal  c^  then  a counter  is  compared  with  the  number  four.  If  four  cases 
have  been  tried  then  an  "unsuccessful  output"  is  produced  stating  this 
fact.  This  kind  of  result  has  not  yet  appeared  in  practice  (April  11, 
1978),  however,  the  possibility  has  been  allowed  for.  If  four  cases 
have  not  yet  been  tried  c-j  is  set  equal  to  c?  and  the  program  returns  to 
that  case  for  another  set  of  iterations. 

Newton-Raphson  Iterations.  The  flowchart  for  this  section  appears  in 
Figure  2.  The  iteration  counter,  t-j , allows  up  to  30  trials.  It  has 
been  the  usual  experience  with  this  program  that  if  convergence  has  not 
been  reached  in  30  iterations,  it  will  not  be  reached  by  going  beyond 
this  point.  Commonly,  convergence  is  attained  with  between  10  and  20 
trials.  It  would  take  less,  but,  for  the  cautiously  slow  start  enforced 
by  the  weighting  factor,  k. 
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Figure  2b 
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Figure  2d 
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Initially,  the  seven  functions  which,  with  slight  variations,  are 
common  to  all  four  cases,  are  evaluated  using  the  previously  determined 
n..  The  variations  required  for  the  different  cases  are  achieved  by 
using  two  constant  multipliers  b and  d,  which  are  set  independently 
equal  to  zero  or  unity  depending  upon  the  case.  A third  constant,  e,  is 
the  subscript  for  A , the  Lagrangian  multiplier  for  the  element  carbon 
in  cases  III  and  IV  when  solid  carbon  is  not  present.  Since  A is 
always  last  in  these  cases,  e is  set  equal  to  n,  the  size  of  the  appropriate 
Jacobian.  Next,  via  an  "on  c^"  statement,  the  appropriate  remaining 
functions  for  cases  I,  III,  and  IV  are  assigned.  This  is  accomplished 
with  the  user-defined  BASIC  functions  fnw  and  fnc. 

Following  the  assignment  of  the  functions,  the  elements  of  the 
upper  left  5x5  submatrix  of  the  Jacobian  are  evaluated.  As  stated 
earlier,  these  depend  upon  whether  the  problem  is  for  constant  pressure 
or  constant  volume. 

Once  evaluated,  the  Jacobian  is  inverted  and  unweighted  corrections 
calculated.  The  next  part  of  the  diagram  introduces  the  weight  factor, 
k,  which  is  applied  only  to  the  five  gases  and  then  up  to  only  trial  15. 
Following  the  correction  of  all  n^ , the  number  of  moles  of  gas  is 
calculated. 

The  next  operation  is  a loop  testing  all  n^  against  the  stopping 
criterion.  If  successful  for  all,  s,  the  counter  of  iteration  sets  is 
incremented  by  one  and  the  calculation  proceeds  to  testing  the  condensed 
phase  case  found  against  the  one  assumed. 
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Case-Consistency  Check.  In  Figure  3,  this  section  begins  by  calculating 
mole  fractions  for  the  five  gases.  These  are  used  to  calculate  equi- 
librium constants  for  the  three  independent  reactions  and  the  partial 
pressure  of  water.  As  stated  earlier  these  are  used  to  check  the  final 
solution  found  by  minimization  of  the  reduced  relative  Planck  function. 

(See  the  subsection  on  Checking  the  Estimates  for  Stopping.)  Although 
only  k^,  for  the  di sproporti onation  reaction,  is  used  in  this  section 
they  are  all  calculated  together. 

The  determination  of  the  case  found  by  the  converged  iterations 
depends  upon  the  case  that  was  assumed.  If  that  was  case  III  or  IV  then 
the  conservation  of  the  element  carbon  was  enforced.  In  these  circum- 
stances, by  comparing  the  equilibrium  constant  for  disproportionation 
(chemical  equation  C)  calculated  with  the  mole  fractions  with  that 
calculated  from  the  Gibbs  energies  of  formation  the  need  for  solid 
carbon  can  be  determined.  On  the  other  hand,  if  carbon  is  assumed 
present  then  the  defined  function  fnc  indicates  whether  its  presence  has 
been  verified. 

For  liquid  water  cases  I and  III  go  together  by  assuming  its  presence. 
For  those  cases  the  sign  of  n^  is  the  indicator.  For  cases  II  and  IV, 
n-j  is  not  calculated,  therefore,  the  partial  pressure  of  water  is  checked 
against  the  vapor  pressure. 

The  value  of  c 2 is  built  up  by  an  initial  assignment  and  addition 
as  the  program  works  its  way  through  this  section.  Upon  reaching  its 
final  value,  c ^ is  compared  with  c-j  as  stated  earlier.  Before  returning 
to  the  appropriate  case,  if  that  is  necessary,  c-j  is  set  equal  to  c 
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Also  it  is  necessary  to  clear  the  eigth  and  ninth  rows  and  columns  of 
the  Jacobian  preparatory  to  inserting  the  appropriate  values  for  the 
next  case. 
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DISCUSSION  OF  THE  LISTING 

Appendix  1 gives  the  list ‘of  symbols  for  the  mathematical  develop- 
ment in  the  text  and  Appendix  2 is  for  the  source  program.  Wherever 
possible,  the  same  variable  symbol  has  been  used  in  the  text  and  in  the 
computer  program.  To  avoid  cluttering  the  text,  variables  used  there 
often  have  fewer  subscripts  than  the  corresponding  symbols  of  the 
program. 

The  listing,  which  is  Appendix  3,  is  given  in  the  full  ASCII 

character  set.  As  stated  earlier,  it  is  dependent  upon  the  BASIC 

* 

compiler  used  by  the  Computer  Sciences  Corporation.  At  this  writing 
(April  11,  1978)  there  is  no  standard  for  the  BASIC  language  as  there  is 
for  FORTRAN  [22].  On  the  other  hand  almost  all  BASIC  compilers  provide 
reasonably  useful  diagnostics  and  internal  editing  commands  so  that 
translation  to  another  compiler  should  not  be  a great  chore. 

Use  of  the  full  ASCII  set  permits  use  of  upper  and  lower  case 
letters.  This  has  been  used  to  differentiate  vectors  and  matrices. 

The  former,  even  when  the  unity  dimension  is  explicitly  indicated,  are 
given  by  lower  case  letters,  as  are  unsubscripted  variables.  Matrices 
are  in  upper  case  as  are  elements  of  matrices. 

Sections  are  separated  by  remark  statements  with  leading  single 
quote  characters.  Some  statements  are  amplified  by  remarks  after  a 


* 

A particular  commerical  source  is  identified  in  this  paper  in  order 
to  explain  the  procedure  adquately.  In  no  case  does  such  identification 
imply  recommendation  or  endorsement  by  the  National  Bureau  of  Standards, 
nor  does  it  imply  that  the  service  is  necessarily  the  best  for  the 
purpose. 
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single  quote  character.  Loops  have  been  doubly  indented  in  order  to 
catch  the  eye,  and  inner,  nested  loops  are  each  indented  again.  The 
listing  is  packed,  i.e.,  blanks  have  been  omitted  in  the  program  state- 
ments except  for  explanatory  matter  and  output  text.  The  experienced 
BASIC  programmer  reads  this  with  ease,  and  a some  storage  space  is 
saved. 
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INFORMATION  FOR  USING  THE  PROGRAM 

Input.  After  signing  on  at  the  terminal  the  program  COLGAS  is  LOADed 
and  RUN.  The  following  queries  will  result. 

The  program  first  asks  whether  the  problem  is  to  be  for  constant 
pressure  or  volume  and  whether  a short  or  long  printout  is  desired.  The 
long  printout  is  generally  only  used  when  attempting  to  rerun  a non- 
converging problem.  The  answer  to  these  joined  questions  is  therefore 
generally  "p,s"  or  "v,s".  For  the  long  printout  "s"  is  replaced  by  "1". 

The  next  question  concerns  the  temperature  and  pressure  (or  volume). 
The  answers,  again  separated  by  a comma,  are  in  kelvins  and  atmospheres 
(or  liters).  Alternatively,  degrees  Celsius  and/or  megapascals  may  be 
used  if  the  values  are  preceded  by  an  asterisk  without  separation  by 
commas  or  blanks. 

The  concentrations  of  the  assumed  gases  are  requested  next,  and  the 
response  is  to  be  in  mole  or  volume  (decimal)  fractions.  The  answers 
are  typed  in  the  order  CH^,  CO,  C02,  H20  with  separating  commas.  If 

any  of  these  are  not  present,  zeros  must  be  entered  in  the  appropriate 
places.  If  the  input  concentrations  do  not  satisfy  certain  requirements, 
the  operator  will  be  requested  to  try  again. 

The  computer  program  does  not  consider  problems  in  which  all  three 
elements,  C,  H,  and  0 are  not  present,  i.e.,  if  CH^  + H2,  H2  + H20,  or 
CO  + C02  are  used  as  starting  mixtures.  In  the  first  two  cases  only 
decomposition  of  the  compounds  can  occur;  in  the  third  only  dispro- 
portionation. (It  is  assumed  that  the  decomposition  of  H2  and  CO  are 

negligible.)  Such  equilibria  are  calculated  by  hand  directly  and 
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without  iteration.  It  is  also  true  that  a mixture  of  CO^  and  F^O  does 
not  undergo  appreciable  reaction. 

At  this  stage,  input  for  the  short  output  is  complete.  If  the  long 
output  was  chosen  the  following  extra  input  is  requested. 

The  opportunity  to  choose  the  starting  case  will  be  next  presented. 
The  response  is  "yes"  or  "no".  If  "yes"  is  typed,  the  (Arabic)  number 
of  the  preferred  case  is  requested. 

Finally,  the  operator  is  to  decide  if  a preferred  set  of  initial 
estimates  is  to  be  used;  "yes"  or  "no".  If  "yes",  the  preferred  values 
are  requested  and  these  are  typed  in  the  same  fashion  as  the  assumed  gas 
concentrations.  Estimates  for  only  the  five  gases  are  to  be  entered. 

This  is  the  end  of  the  operator  input. 

Output.  All  input  data  are  reprinted  as  an  echo  check.  After  the 
input  gas  concentrations  have  been  echoed,  the  unchosen  quantity,  volume 
or  pressure,  is  given.  Because  the  number  of  moles  of  gas  almost 
always  changes  in  the  final  equilibrium,  all  extensive  quantities  are 
quoted  on  the  basis  of  the  number  of  moles  of  the  assumed  gas.  Although 
pressure  is  not  an  extensive  quantity  it  too  can  only  be  quoted  on  this 
basis  when  a constant  volume  problem  is  calculated  because  the  number  of 
moles  of  gas  in  the  input  volume  has  not  been  requested. 

A long  printout  gives  at  this  point  the  assumed  case  and  the  results 
of  each  iteration  until  the  program  terminates  successfully  or  unsuccess- 
fully. In  the  latter  situation  a statement  about  unsuccessful  convergence 
is  followed  by  the  last  set  of  estimates.  In  any  case  the  mole  fractions 
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are  next  and  these  are  followed  by  the  number  of  moles  of  gas  and  the 
new  volume  (or  pressure). 

If  liquid  water  or  solid  carbon  are  produced  a pertinent  statement 
is  made  at  this  point.  If  the  temperature  is  above  the  critical  value 
for  water,  such  a statement  is  made  next. 

Normally,  that  is  the  end  of  the  printout.  However,  if  the  run  was 
unsuccessful  the  independent  equilibria  will  not  have  been  satisfied  and 
statements  to  this  effect  will  be  the  last  output. 
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THE  PHASE  RULE 

For  the  restricted  sets  of  problems  considered  here,  the  phase  rule 
is  simple  to  apply.  It  is  given  by  equation  (48). 

F = C - P + 2 (48) 

The  number  of  degrees  of  freedom,  F,  is  the  number  of  independent 
choices  of  variables  possible  with  a given  system  which  has  P phases  and 
C components.  In  order  to  use  the  traditional  symbols  for  the  phase 
rule,  P has  been  used  here,  although  it  has  been  used  before.  Because 
this  applies  only  to  this  section,  confusion  should  not  result.  There 
are  three  possible  phases,  gas,  liquid  water,  and  solid  carbon.  The 
number  of  components  in  these  problems  is  equal  to  the  number  of  chemical 
elements,  three.  This  simplifies  equation  (48). 

F = 5 - P (49) 

In  case  IV  (see  the  section  on  General  Statement  of  the  Algorithm) 
there  is  only  the  gas  phase,  and  there  must  be  four  variables  available 
for  choice.  Since  two  are  assumed  to  be  temperature  and  pressure  (or 
volume),  there  are  two  left  for  concentrations . Thus,  if  a triangular 
diagram  is  used  for  the  concentrations  of  the  elements  carbon,  hydrogen, 
and  oxygen,  the  possible  elemental  concentrations  would  be  represented 
by  an  area.  See  Figure  4.  If  two  atomic  compositions  are  chosen,  the 
third  is  fixed  by  the  requirement  that  the  percentages  add  to  100. 

On  the  other  hand  for  an  equilibrium  that  has  either  solid  carbon 
or  liquid  water  (but  not  both)  present,  there  is  only  one  degree  of 


Figure 
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freedom,  and  this  must  be  represented  by  a line.  Choosing  the  percentage 
for  just  one  element  fixes  both  of  the  others  if  one  condensed  phase  is 
present.  Also,  because  the  condensed  phase  always  has  a fixed  com- 
position, the  variable  composition  refers  to  the  gas  phase.  Referring 
to  Figure  3,  the  curved  line  shows  the  gas  compositions  in  equilibrium 
with  solid  carbon  at  700  kelvins  and  1 atmosphere.  The  curve  has  been 
drawn  by  connecting  points  calculated  with  COLGAS  which  turned  out  to  be 
case  II.  If  a gas  mixture  contains  15  atomic  percent  carbon  and  is  in 
equilibrium  with  solid  carbon  the  the  gas  could  contain  also  about  48 
atomic  percent  hydrogen  and  37  percent  oxygen.  Alternatively,  it  could 
contain  81  percent  hydrogen  and  4 percent  oxygen.  The  two  possibilities 
are  at  the  intersections  of  the  15  percent  carbon  line  with  the  calculated 
curve. 

A valuable  use  of  the  phase  diagram  in  Figure  4 is  to  show  potential 
carbon  precipitation  for  an  assumed  test  gas  mixture.  One  needs  merely 
to  calculate  the  elemental  composition  from  the  compound  composition  and 
locate  that  point  on  the  diagram.  If  it  falls  in  the  area  above  the 
curve,  carbon  will  precipitate.  Furthermore,  extending  the  line  through 
the  apex  (for  c = 100%)  and  the  point  to  the  curve  will  give  the  elemental, 
but,  not  the  compound  composition  of  the  gas  in  equilibrium  with  the 
carbon.  It  should  be  noted  that  gas  mixtures  made  from  the  five  postulated 
substances  will  have  elemental  compositions  within  the  region  bounded  by 
connecting  the  compositions  of  all  of  the  substances  with  straight 
lines.  The  lines  should  connect  CH^,  CO,  CO2,  F^O,  F^,  and  return  to 
CH^  in  that  order. 
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Figure  4 represented  a situation  above  the  critical  temperature  of 
water.  Figure  5 illustrates  the  situation  for  500K  and  50  atmospheres. 
The  diagram  is  rather  more  complicated.  Although  not  always  indicated, 
a gas  phase  is  present  for  all  regions  of  the  diagram.  There  are  two 
invariant  points  and  two  regions  within  which  gas  compositions  will 
change  to  one  or  the  other  resulting  in  condensed  water  and  carbon. 
There  are  two  regions  within  which  the  change  will  result  in  carbon 
only,  and  two  more  from  which  only  water  will  condense.  Note  the  two 
small  regions  within  which  no  condensation  will  occur,  although  five 
gaseous  substances  may  change  their  relative  amounts. 
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EXAMPLES 

Examples  A through  D are  representative  of  cases  I through  IV  and 
are  successful  short  printout  results.  Example  B,  for  700K  and  1 
atmosphere,  corresponds  to  the  conditions  for  the  phase  diagram  of 
Figure  3.  The  original  elemental  proportions  are  .6:1. 6:. 8 and  the 
atomic  fractions  are  .2,  .533,  and  .267.  The  line  extended  from  the 
apex  through  this  point  on  the  diagram  intersects  the  curve  for  solid 
carbon  at  approximately  .13,  .58,  and  .29.  The  equilibrium  mole  fractions 
for  the  five  substances  give  for  the  atomic  fractions  in  the  order  C, 

H,  0:  .1288,  .5809,  and  .2904. 

Example  E is  for  a problem  which  did  not  converge  during  the  short- 
printout  run.  The  last  set  of  estimates  suggests  that  the  difficulty 
involved  the  amounts  of  CO  and  CO^  which  appeared  to  be  continually 
going  negative.  Often  in  a situation  of  this  kind  it  was  found  that  the 
initial  estimates  for  such  substances  were  not  small  enough.  Consequently, 
for  the  long  printout  run,  example  F,  initial  estimates  were  chosen 
equal  to  the  last  set  from  the  short  run  except  for  the  CO  and  CO2  which 


were  made  smal ler. 
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Example  A 

col  gas  13:14  04/07/78 

Is  this  a constant  pressure  (p),  or  a constant  volume  (v)  problem,  and 
do  you  want  a short  (s)  or  a long  (1)  printout?p,s 

What  are  the  temperature  (kelvins)  and  pressure  (atmospheres)?500 , 1 00 
Temperature  = 500  K = 226.85  C 
Pressure  = 100  atm  = 10.1325  MPa 

What  are  the  assumed  concentrations  in  mole  fractions  of  the  gases  in 
the  order:  CH^,  CO,  CO^,  H^,  and  H^O? . 1 , . 1 , . 1 , . 4 , . 3 
.1  .1  .1  .4  .3 

The  initial  volume  = .410284  liters  per  mole  of  original  gas 

Gas  contains  C,  H,  and  0 atoms  in  the  proportions  .3  , 1.8  , and  .6 

Equilibrium  mole  fractions  in  the  order  CH^,  CO,  CO2,  H^O  are: 
.7301  .0000  .0077  .0016  .2606 


The  number  of  moles  of  gas  = .2074  per  mole  of  original  gas 

The  volume  = 8.50928e~2  liters  per  mole  of  original  gas 

Solid  carbon  is  present  in  the  amount  of  .1470  moles  per  mole  of  original 

gas 

Liquid  water  is  present  in  the  amount  of  .5428  mole  per  moles  of  original 
gas 

now  at  end 
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Example  B 

colgas  13:16  04/07/78 

Is  this  a constant  pressure  (p),  or  a constant  volume  (v)  problem,  and 
do  you  want  a short  (s)  or  a long  (1)  printout?p,s 

What  are  the  temperature  (kelvins)  and  pressure  (atmospheres)?700 , 1 

Temperature  = 700  K = 426.85  C 
Pressure  = 1 atm  = .101325  MPa 

What  are  the  assumed  concentrations  in  mole  fractions  of  the  gases  in 
the  order:  CH4,  CO,  C02,  H2,  and  H20?. 2, . 2, . 2, . 2, . 2 
.2  .2  .2  .2  .2 

The  initial  volume  = 57.4398  liters  per  mole  of  original  gas 

Gas  contains  C,  H,  and  0 atoms  in  the  proportions  .6  , 1.6  , and  .8 

Egui librium  mole  fractions  in  the  order  CH^,  CO,  C02,  H^,  H20  are:  .1693 
.0076  .2344  .1376  .4511 

The  number  of  moles  of  gas  = .8626  per  mole  of  original  gas 

The  volume  = 49.5486  liters  per  mole  of  original  gas 

Solid  carbon  is  present  in  the  amount  of  .2452  moles  per  mole  of 
original  gas 

Above  the  critical  temperature  for  water 
now  at  end 


58 


Example  C 

col  gas  13:18  04/07/78 

Is  this  a constant  pressure  (p),  or  a constant  volume  (v)  problem,  and 
do  you  want  a short  (s)  or  a long  (1)  printout?p,s 

What  are  the  temperature  (kelvins)  and  pressure  (atmospheres)?500 ,50 
Temperature  = 500  K = 226.85  C 
Pressure  = 50  atm  = 5.06625  MPa 

What  are  the  assumed  concentrations  in  mole  fractions  of  the  gases  in 
the  order:  CH^,  CO,  C02>  H^,  and  h^O?. 05, . 1 , . 05, . 5, . 3 
.05  .1  .05  .5  .3 

The  initial  volume  = .820569  liters  per  mole  of  original  gas 

Gas  contains  C,  H,  and  0 atoms  in  the  proportions  .2  , 1.8  , and  .5 

Equilibrium  mole  fractions  in  the  order  CH^,  CO,  CO^,  are:  .4730 

.0000  .0012  .0047  .5211 

The  number  of  moles  of  gas  = .4218  per  mole  of  original  gas 
The  volume  = .34611  liters  per  mole  of  original  gas 

Liquid  water  is  present  in  the  amount  of  .2792  moles  per  mole  of  original 
gas 

now  at  end 
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Example  D 

col  gas  14:22  05/03/78 

Is  this  a constant  pressure  (p),  or  a constant  volume  (v)  problem,  and 
do  you  want  a short  (s)  or  a long  (1)  printout?v,s 

What  are  the  temperature  (kelvins)  and  volume  ( 1 i ters )? 1255,1.514 
Temperature  = 1255  K = 981.85  C 
Volume  = 1.514  1 iters 

What  are  the  assumed  concentrations  in  mole  fractions  of  the  gases  in 
the  order:  CH^,  CO,  CO2,  , and  H^O?. 05, . 18, . 12, . 25, . 4 
.05  .18  .12  .25  .4 

The  initial  pressure  - 68.0194  atm  per  mole  of  original  gas 

Gas  contains  C,  H,  and  0 atoms  in  the  proportions  .35  , 1.5  , and  .82 

Equilibrium  mole  fractions  in  the  order  CH^,  CO,  C02#  E^O  are: 
.0169  .1902  .1218  .3342  .3368 

The  number  of  moles  of  gas  = 1.0640  per  mole  of  original  gas 
The  pressure  = 72.37  atm  = 7.33289  MPa  per  mole  of  original  gas 

Above  the  critical  temperature  for  water 
now  at  end 
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Example  E 

colgas  14:14  05/03/78 

Is  this  a constant  pressure  (p),  or  a constant  volume  (v)  problem,  and 
do  you  want  a short  (s)  or  a long  (1)  printout?p,s 

What  are  the  temperature  (kelvins)  and  pressure  (atmospheres)?500,50 
Temperature  = 500  K = 226.85  C 
Pressure  = 50  atm  = 5.06625  MPa 

What  are  the  assumed  concentrations  in  mole  fractions  of  the  gases  in 
the  order:  CH4,  CO,  C02,  H2,  and  H20?. 12,0, .01 ,.75,. 12 
.12  0 .01  .75  .12 

The  initial  volume  = .820569  liters  per  mole  of  original  gas 

Gas  contains  C,  H,  and  0 atoms  in  the  proportions  .13  , 2.22  , and  .14 
30  iterations  did  not  converge! 

The  results  at  this  stage  are  as  follows 
Case  IV 

The  last  set  of  estimates  is,  moles  and  multipliers: 

.13  l.e-9  l.e-9  .71  .141.79487  54.3204  13.1832 

Equilibrium  mole  fractions  in  the  order  CH^,  CO,  C02>  H2>  H20  are:  .1327 
.0000  .0000  .7245  .1429 

The  number  of  moles  of  gas  = .9800  per  mole  of  original  gas 

The  volume  = .804157  liters  per  mole  of  original  gas 

Equilibrium  constants  do  not  check  out 
First  equilibrium  constant  = 7.24886e~3  vs.  .197184 
Second  equilibrium  constant  = 1.20341e+10  vs.  19534.9 
Third  equilibrium  constant  = 1.61 777e- 9 vs.  5.10204e~8 


now  at  end 
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Example  F 

col  gas  14:16  05/03/78 

Is  this  a constant  pressure  (p),  or  a constant  volume  (v)  problem,  and 
do  you  want  a short  (s)  or  a long  (1)  printout?p,l 

What  are  the  temperature  (kelvins)  and  pressure  (atmospheres)?500,50 
Temperature  = 500  K = 226.85  C 
Pressure  = 50  atm  = 5.06625  MPa 

What  are  the  assumed  concentrations  in  mole  fractions  of  the  gases  in 
the  order:  CH^,  CO,  CO2,  H^,  and  F^O?. 12,0, .01 ,.75,. 12 
.12  0 .01  .75  .12 

The  initial  volume  = .820569  liters  per  mole  of  original  gas 

Gas  contains  C,  H,  and  0 atoms  in  the  proportions  .13  , 2.22  , and  .14 
Do  you  want  to  choose  the  starting  case?no 
no 

Do  you  want  to  provide  your  own  initial  estimates?yes 
yes 

What  are  they  in  the  same  order  as  above?  .13  l.e-15  l.e-14  .71 

.14 

0 .13  l.e-15  l.e-14  .71  .14 

0 0 0 
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Case  IV 


1 

.13 

1 . 00465e- 1 5 

1.01427e-14 

.71 

. 14 

-1. 

79487 

54.3204 

13. 1832 

2 

. 13 

1 . 023e- 1 5 

1.071053-14 

.71 

. 14 

-1. 

79487 

54.3204 

13. 1832 

3 

.13 

1 . 0628e- 1 5 

1. 1989e-14 

.71 

.14 

-1. 

79487 

54.3204 

13. 1832 

4 

.13 

1 . 1 2872e~ 1 5 

1 .4278e-14 

.71 

. 14 

-1. 

79487 

54.3204 

13. 1832 

5 

.13 

1 . 21 986e~ 1 5 

1.78602e-14 

.71 

. 14 

-1. 

79487 

54.3204 

13. 1832 

6 

.13 

1.32744e-15 

2. 28501 e~ 14 

.71 

. 14 

-1. 

79487 

54.3204 

13. 1832 

7 

.13 

1 .43493e-15 

2.89084e-14 

.71 

. 14 

-1. 

79487 

54.3204 

13.1832 

8 

.13 

1.52342e-15 

3. 5031 3e- 1 4 

.71 

.14 

-1. 

79487 

54.3204 

13. 1832 

9 

.13 

1 . 581 22e- 1 5 

3. 98772e- 14 

.71 

. 14 

-1. 

79487 

54.3204 

13. 1832 

10 

. 13 

1 . 60992e- 1 5 

4.26874e-14 

.71 

. 14 

-1. 

79487 

54.3204 

13. 1832 

11 

.13 

1 . 62027e- 1 5 

4.38108e-14 

.71 

.14 

-1. 

79487 

54.3204 

13. 1832 

12 

.13 

1 . 62284e- 1 5 

4.41047e-14 

.71 

. 14 

-1. 

79487 

54.3204 

13. 1832 

13 

.13 

1 . 62325e- 1 5 

4.41 521 e- 1 4 

.71 

. 14 

-1. 

79487 

54.3204 

13. 1832 

14 

.13 

1 . 62328e- 1 5 

4.41 564e~14 

.71 

. 14 

-1. 

79487 

54.3204 

13. 1832 

15 

.13 

1 . 62329e- 1 5 

4. 41 565e- 14 

.71 

.14 

-1. 

79487 

54.3204 

13. 1832 

Equi librium 

mole  fractions 

in  the  order  CH^, 

CO, 

co2, 

h2. 

H20  are: 

.132653  1 . 65641 e- 1 5 4.50576e-14  .72449  .142857 

The  number  of  moles  of  gas  = 1.  per  mole  of  original  gas 
The  volume  = .820569  liters  per  mole  of  original  gas 
now  at  end 
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SYMBOLS  - TEXT 

Helholtz  energy  of  the  system 

Number  of  components 

Number  of  gram  atoms  of  element  j 

Number  of  degrees  of  freedom 

Function  actually  minimized  for  constant  volume  problems 
Function  actually  minimized  for  constant  pressure  problems 

Gibbs  energy  of  the  system 

Gibbs  energy  of  the  chemical  elements 

Jacobian  of  the  functions  set  equal  to  zero  for  minimization 
Inverse  of  J 

Function  expressing  conservation  of  element  j 
Number  of  chemical  species 

Number  of  phases  (discussion  of  phase  rule  only) 

Pressure  of  the  system 

Number  of  independent  chemical  reactions  (subsection  Checking  the 
Estimates  for  Stopping  Only) 

Universal  gas  constant 
Thermodynamic  temperature 
Volume  of  the  system 

Coefficients  in  the  empirical  equation  for  In  k- 

Coefficients  in  the  empirical  equation  for  the  natural  logarithm 

of  the  vapor  pressure  of  water 

Number  of  gram  atoms  of  carbon  in  the  system 

Number  of  chemical  elements  in  the  system 

Value  at  the  t-th  iteration  of  the  functions  set  equal  to  zero 

for  minimization 

Number  of  gaseous  substances 

Number  of  gram  atoms  of  hydrogen  in  the  system 

Equilibrium  constant  of  formation  of  A B 

Equilibrium  constant  of  formation  of  tne^i-th  substance 

Value  of  weighting  factor  k for  t-th  trial 

Natural  logarithm 

Maximum  exponent  of  temperature  in  equations  for  equilibrium 

constants  of  formation  and  vapor  pressure  of  water 

Number  of  moles  of  solid  carbon 

Total  number  of  moles  of  gas 

Number  of  moles  of  the  i-th  substance 

Number  of  gram  atoms  of  oxygen  in  the  system 

Partial  pressure  of  substance  i 

Saturated  vapor  pressure  of  water 

Total  number  of  substances 
Number  of  trial  or  interation 

Coefficients  and  subscripts  in  chemical  equation  of  formation 
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AG° 

rA 

ag 


A B 

x y 


with  all 


Vector  representing  the  variables  being  estimated  at 
iteration  number  t 

Gibbs  energy  of  formation  of  one  mole  of 
substances  in  the  standard  states 
Function  of  A which  is  minimized 
Function  of  G which  is  minimized 
Lagrangian  multiplier 

Chemical  potential,  partial  molar  Gibbs  energy  of  the 
i-th  substance 


65 


I ( i , j ) 
J(i , j) 
N(i , tl ) 

a$ 

a,b,d,e 

b$ 

c 

c$ 

c(i,D 

c(i  )$ 

cl 

c2 

d$ 

dl 


f(i,D 

fnc 

fnw 

9 

g(ti ) 

h 

k 


k k k 
4*5*6 


Hi) 


n 


P 

p$ 

p( tl ) 
pO 


t 

t$ 

tO , t9 


APPENDIX  2 

SYMBOLS  - PROGRAM 

Inverse  of  the  Jacobian, J 
Jacobian 

Number  of  moles  of  the  i-th  substance  for  the 
t-th  trial 

"Yes"  or  "no"  to  query  about  choosing  case 
Constants  used  to  differentiate  functions  for  the 
cases 

"Yes"  or  "no"  to  query  about  choosing  initial  estimates 
Number  of  gram  atoms  of  carbon  in  the  system 
"p"  or  "v"  for  constant  pressure  or  constant  volume 
problem 

Unweighted  correction  to  the  i-th  variable 

Case  number  i as  a string  variable 

Case  number  during  an  iteration  set 

Case  number  found  as  a result  of  iterations 

"1"  or  "s"  for  length  of  printout  desired 

Relative  difference  between  saturated  vapor  pressure 

of  water  and  partial  pressure  of  water  calculated  from 

mole  fraction  found  by  iteration 

The  i-th  function  being  set  equal  to  zero 

Carbon  function 

Water  function 

Original  number  of  moles  of  gas 
Number  of  moles  of  gas 

Number  of  gram  atoms  of  hydrogen  in  the  system 
Weighting  factor  for  applying  corrections 
Equilibrium  constants  for  the  three  independent 
reactions  calculated  from  the  logarithms  of  the 
equilibrium  constants  of  formation 
Equilibrium  constants  for  the  three  independent 
reactions  calculated  from  mole  fractions  found 
by  iteration 

Natural  logarithm  of  equilibrium  constant  of 

formation  of  the  i-th  substance 

Size  of  the  pertinent  Jacobian 

Number  of  gram  atoms  of  oxygen  in  the  system 

Pressure  in  atmospheres 

Pressure  as  a string  variable 

Pressure  of  a constant  volume  system  at  the  tl-th 
trial 

Initial  pressure  in  a constant  volume  system 

Universal  gas  constant 

Number  of  sets  of  iterations  performed 

Temperature  in  kelvins 

Temperature  as  a string  variable 

Starting  and  ending  indexes  for  the  iteration 

loop 


cm  ro  <3- 
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t-.  Number  of  an  iteration 

* Square  of  t 

Cube  of  t 
Fourth  power  of  t 
Volume 

vO  Initial  volume  in  a constant  pressure  system 

x(i)  Mole  fraction  of  the  i-th  gaseous  substance 
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APPENDIX  3 
LISTING 

TOO'  colgas  WSHorton  9/7/77 

110'  Calculates  equilibrium  composition  for  the  CH4,  CO,  C02,  H2,  H20 
system 

120'  by  minimizing  the  Gibbs  energy  for  a constant  pressure  problem  or  the 
130'  Helmolz  energy  for  a constant  volume  problem. 

140' 

150'  Initial  Statements 
160' 

170width80 

180randomize 

190dimg( 150) ,1(15,15) ,J(15,15) ,N(15,150) ,p(150) 

200deffnw=2*n(6,tl-l )+n(7,tl-l)-l(6) 

21 Odef fnc=c-n( 1 ,tl-l)-n(2,tl-l)-n(3,tl-l) 

220r=. 08205686 ' V-zero  = 24. 41383e-3mt3molet-l 
230c ( 1 )$=" I" 

240c(2)$="II" 

250c(3)$=" III" 

260c(4)$=" IV" 

270' 

280'  Input 
290' 

300print"Is  this  a constant  pressure  (p),  or  a constant  volume  (v)  problem," 
310print"and  do  you  want  a short  (s)  or  a long  (1)  printout"; 

320 i nputc$ ,d$ 

330print"What  are  the  temperature  (kelvins)  and"; 

340 i fc$="p"orc$="P"then380 
350print"  volume  (liters)"; 

360 i nputt$ , v 
370goto400 

380print"  pressure  (atmospheres)"; 

390 i nputt$ ,p$ 

400 i f t$( 1 , 1 )="*" then430 
41 0t=val (t$) 

420goto440 

430t=val (t$ (2)) +273. 15 

440pri nt"Temperature  ="t;"K  =" ; t- 273 . 15;"C" 

450 i fc$="p"orc$="P"then  480 
460pri nt"Vol ume  ="v;"l iters" 

470goto530 

480i fp$( 1 ,l)="*"then510 
490p=val (p$) 

500goto520 

510p=val(p$(2))/. 10325 

520pri nt"Pressure  ="p;"atm  =" ;p*0. 101325;"MPa" 

530pri nt 

540pri nt"What  are  the  assumed  concentrations  in  mole  fractions  of  the  gases" 
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550print"in  the  order:  CH4,  CO,  C02,  H2,  and  H20" ; 

560 i nputn( l,0),n(2,0),n(3,0),n(4,0),n(5,0) 

570g=g(0)=n( 1 , 0)+n(2 ,0)+n(3 ,0)+n(4,0)+n(5 ,0) 

580c=n(l ,0)+n(2,0)+n(3,0) 

590h=4*n(l ,0)+2*n(4,0)+2*n(5,0) 

600o=n(2 ,0)+2*n(3 ,0)+n(5 ,0) 

61  Oi  fc<t<>0andh<t<>0ando<t<>0then650 

620print"It  is  necessary  to  have  compounds  containing  carbon,  oxygen," 
630print"and  hydrogen.  Please  try  again." 

640goto540 

650 i f n ( 1 , 0)+n(2 ,0)+n(4 ,0)<t<>0then680 

660pri nt"Carbon  dioxide  and  water  alone  do  not  react  significantly" 

670goto3850'  print, end 

680fori=l to5 

690  printn(i,0); 

700  i fn( i ,0)<>0then720 

710  n(i,0)=le-7 

720nexti 

730pri nt 

740pri nt 

750 i fc$="p"orc$="P"then810 
760a=0 

770p=p0=r*t/v 

780p(0)=g(0)*p0 ' Initial  estimate  of  pressure 

790print"The  initial  pressure  ="p(0);"atm  per  mole  of  original  gas" 

800goto840 

810a=l 

820v0=r*t/p'  Initial  volume 

830print"The  initial  volume  ="v0;"l iters  per  mole  of  original  gas" 

840pri nt 

850print"Gas  contains  C,  H,  and  0 atoms  in  the  proportions'^;" ,"h;" , and"o 
860ifd$="s"ord$="S"then970 

870print"Do  you  want  to  choose  the  starting  case"; 

880 i nputa$ 

890pri nta$ 

900 i fa$="no"then960 

91 Opri nt"What  case  number"; 

920 i nputcl 

930pri nt"Case  number"cl 
940pri nt 
950goto990 
960pri nt 

970c 1 =0 1 Causes  bypass  of  line  1530 
980 i fd$="s"ord$="S"thenl 180 

990print"Do  you  want  to  provide  your  own  initial  estimates"; 

1 000 i nputb$ 

1 01 Opri ntb$ 

1 020pri nt 

1 030i fb$="yes" then  1 050 
1 040gotol 090 
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1050pri nt"What  are  they  in  the  same  order  as  above"; 

1 060 i nputn(l ,0),n(2,0),n(3,0),n(4,0),n(5,0) 

1070printn(l ,0);n(2,0);n(3,0);n(4,Q);n(5,0) 

1080print 
1090pri nt"  0"; 

1 1 00g(0)=0 
1 1 10fori=l to5 
1120  printn(i,0); 

1130  ifn(i ,0)<>0thenl 150 
1140  n(i ,tl)=10i(-7-int(4*rnd)) 

1150  g(0)=g(0)+n(i ,0) 

1 160nexti 
1 170print 

1 180n(6,0)=n(7,0)=n(8,0)=n(9 ,0)=0 
1190' 

1200'  Thermodynamic  Data 
1210' 

1220t2=t*t 
1 230t3=t*t2 
1 240t4=t2*t2 

1 250£( 1 )=821 2. 35/ t- . 81 6063*1 og(t)~ 1 . 00801 e~2*t+7. 73804e~6*t2 
1260£(l)=£(l)-3. 06981 e-9*t3+5. 08552e-13*t4 
1 270£(2)=1 3632. 3/t+l . 8055*1  og(t)- . 0023936* t+3.  72836e~7*t2 
1280£(3)=47326.7/t+.  0777872*1  og(t)-2. 80706e-4*t+3. 48747e-8*t2 
1 290£(4)=0 

1 300£(5)=28786. 6/t-. 69731 I*log(t)-.00145957*t+8.3769e-7*t2-1.80285e-10*t3 

1 31 0£(7)=55. 5049-9430. 74/t- . 1 59684*t+3. 07499e-4*t2-3. 0549e-7*t3+l . 24325e~ 1 0*t4 

1320£(6)=1(5)-1(7) 

1330kl=exp(£(2)+£(5)-£(3)-£(4)) 

1340k2=exp(£(l)+£(5)-£(2)-3*£(4)) 

1 350k3=exp(2*£(2)-£(3) ) 

1 360fori=l to5 

1370  m(i )=log(p)-£(i ) 

1 380nexti 
1390' 

1400'  Assign  max.  size  and  constant  elements  of  basic  Jacobian 
1410' 

1420n=9 

1430matJ=zer(n , n) 

1 440 J (2,7 )=J (7,2 )=J (5,7 )=J ( 7 , 5 )=1 

1450J(3 , 7)=J(7 ,3)=J(4 ,6)=J(6 ,4)=J(5 ,6)=J(6 ,5)=2 

1 460 J ( 1 ,6)=J(6,1)=4 

1470' 

1480'  Assign  constants  for  first  set  of  iterations 
1490' 

1 500s=0 
1 51 0t0=l 
1 520t9=30 

1 530onc 1 goto 1 630 , 1 790 , 1 920 , 2080 
1540' 
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1550'  Test  for  possibility  of  liquid  water" 
1560' 

1 5 7 0 i ft<647. 3andp>exp(£(7)) then 1600 

1 580cl=2 

1 590gotol 790 

1 600cl=4 

1610goto2080 

1620‘ 

1630'  Case  I 
1640' 

1 650 i f s>0thenl 680 
1 660 i fd$=" s"ord$="S" then 1 680 
1670pri ntn(6 ,0) ;n(7,0);n(8,0) 

1 680b=e=0 
1 690d=l 
1 700n=8 

1710J(6,8)=J(8,6)=2 
1 720J(7 ,8)=J(8,7)=1 
1 730i fd$="s"ord$="S"thenl 770 
1 740pri nt 
1 750pri nt"Case  I" 

1 760pri nt 

1 770goto2220 1 iterate 
1780' 

1790'  Case  II 
1800' 

1810i fs>0thenl840 
1 820 i fd$="s"ord$="S" then 1840 
1830printn(6,0) ;n(7,0) 

1840b=d=e=0 

1850n=7 

1 860  i fd$="s"ord$="S" then 1 900 
1870pri nt 

1 880pri nt"Case  II" 

1890pri nt 

1900goto2220'  iterate 
1910' 

1920'  Case  III 
1930' 

1 940 i f s>0thenl 970 

1 950 i fd$="s"ord$="S" then 1 970 

1 9 6 0 p r i ntn(6 ,0) ;n(7,0);n(8,0);n(9,0) 

1 970b=d=l 
1 980e=n=9 

1 990J(6 ,8)=J(8,6)=2 
2000J(7 ,8)=J(8,7)=1 

201 0J( 1 ,9)=J(9,1)=J(2,9)=J(9,2)=J(3,9)=J(9,3)=1 
2020 i fd$=" s"ord$="S" then2060 
2030pri nt 

2040pri nt"Case  III" 
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2 0 5 0 p r i nt 

2060goto2220'  iterate 
2070' 

2080'  Case  IV 
2090' 

2100ifs>0then2130 

21 1 Oi fd$="s"ord$="S"then21 30 

2120pri ntn(6 ,0) ;n(7,0);n(8,0) 

21 30b=l 
21 40d=0 
2150e=n=8 

2160J(1 ,8)=J(8,1)=J(2,8)=J(8,2)=J(3,8)=J(8,3)=1 
21 70i fd$="s"ord$="S"then2220 
2180pri nt 

2190pri nt"Case  IV" 

2200pri nt 
2210“ 

2220'  Newton-Raphson  Iterations 
2230' 

2240fortl=t0tot9 

2250  f(l ,l)=log(n(l ,tl-l))-a*log(g(tl-l))+4*n(6,tl-l)+b*n(e,tl-l)+m(l ) 
2260  f(2,l)=log(n(2,tl-l))-a*log(g(tl-l))+n(7,tl-l)+b*n(e,tl-l)+m(2) 

2270  f(3,l)=log(n(3,tl-l))-a*log(g(tl-l))+2*n(7,tl-l)+b*n(e,tl-l)+m(3) 

2280  f(4,l)=log(n(4,tl-l))-a*log(g(tl-l))+2*n(6,tl-l)+m(4) 

2290  f(5,l)=log(n(5,tl-l))-a*log(g(tl-l))+2*n(6,tl-l)+n(7,tl-l)+m(5) 

2300  f (6 , 1 )=4*n( 1 ,tl-l)+2*n(4,tl-l)+2*n(5,tl-l)+d*2*n(8,tl-l)-h 

2310  f(7, 1 )=n( 2 , 1 1 - 1 )+2*n(3,tl-l )+n(5,tl-l )+d*n(8,tl-l  )-o 

2320  onclgoto2330, 2390, 2350, 2380 

2330  f(8,l)=fnw 

2340  goto2390 

2350  f (8 , 1 )=f nw 

2360  f(9,l)=-fnc 

2370  goto2390 

2380  f (8 , 1 )=-fnc 

2390  i fc$="p“orc$="P"then2440 

2400  fori=lto5 

2410  J(i,i)=l/n(i,tl-l) 

2420  nexti 

2430  goto2520 

2440  fori=lto5 

2450  forj=lto5 

2460  i f i=jthen2490 

2470  J(i,J)=-l/g(tl-l) 

2480  goto2500 

2490  J(i , i )=l/n(i , tl-1 )-l/g(tl-l ) 

2500  nextj 

2510  nexti 

2520  matI(n,n)=inv(J(n,n)) 

2530  matc(n,l )=I(n,n)*f(n,l) 

2540  ifd$="s"ord$="S"then2560 
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2550  printtl; 

2560  fori=lton 

2570  ifi>5then261Q 

2580  iftl>15then2610 

2590  k=sin((t.l-8)*&pi/16)/2+l/2 

2600  got.o2620 

2610  k=l 

2620  n(i ,tl)=n(i ,tl-l )-k*c(i , 1 ) 

2630  ifn(i ,tl )>0then2660 

2640  i f i>5then2660 

2650  n(i ,tl )=10t(-7-int(4*rnd)) 

2660  i fd$="s"ord$="S"then2680 

2670  printn(i ,tl ) ; 

2680  nexti 

2690  ifd$="s"ord$="S"then2710 
2700  print 

2710  g(tl)=n(l ,tl)+n(2,tl)+n(3,tl)+n(4,tl)+n(5,tl) 

2720  i fc$="p"orc$="Pllthen2740 
2730  p=p(tl)=g(tl)*pO 
2740  fori=lton 

2750  ifabs(c(i ,1 )/n(i , t 1 ))>5e~5then2780 

2760  nexti 

2770  goto2880 

2780nexttl 

2790tl=tl-l 

2800pri nttl ; " iterations  did  not  converge!" 

2810fori=l to5 

2820  x(i )=n(i , tl )/g(tl ) 

2830nexti 

2840k4=x(2)*x(5)/(x(3)*x(4))'  C02  + H2  = CO  + H20 

2850k5=x( 1 )*x(5)/(x(2)*x(4)i3*p*p) 1 CO  + 3H2  - CH4  + H20 
2860k6=x(2)*x(2)*p/x(3) 1 C02  + C(s)  = 2C0 

2870goto3270 
2880s=s+l 


2890' 

2900'  Case-consistency  checking 
2910“ 

2920fori=l to5 

2930  x(i )=n(i ,tl )/g(tl ) 

2940nexti 

2950k4=x(2)*x(5)/(x(3)*x(4))'  C02  + H2  = CO  + H20 

2960k5-x(l)*x(5)/(x(2)*x(4)i3*p*p)'  CO  + 3H2  = CH4  + H20 
2970k6=x(2)*x(2)*p/x(3) 1 C02  + C(s)  = 2C0 

2980d]=(x(5)*p~exp( 1 (7)))/exp(l (7)) ' H20:  (p.p.  - v.p.)/v.p. 
2990 i fcl>2then3030 
3000 i f f nc>0then3050 


301 0c2=2 

3020goto3060 

3030ifk6>k3then3050 


3040goto301 0 
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3050c2=0 

3060onclgoto3070, 3090, 3070, 3090 
30701 f n(8 , tl )>0then31 20 
3080goto3100 

30901 fx(5)*p>exp(£(7))then3120 
3100c2=c2+2 
31 1 0goto31 30 
31 20c2=c2+l 

31 301 fc2=cl then3340 ' successful  exit 

3 1 4 0 i fs>3then3250 ' unsuccessful  exit;  four  inconsistent  sets  of  iterations 

3150cl=c2 

3160fori=8ton 

3170  forj=lton 

3180  J(i,J)=J(J,i)=0 

3190  nextj 

3200nexti 

3210onclgotol630, 1790, 1920,2080 
3220' 

3230'  Output 
3240' 

3250pri nt"Coul d not  determine  correct  case  with  four  sets  of  iterations" 

3260pri nt 

3270pri nt"The  results  at  this  stage  are  as  follows" 

3280pri nt"Case  "c(cl)$ 

3290pri nt"The  last  set  of  estimates  is,  moles  and  multipliers:" 

3300fori=l ton 
3310  pri ntn( i , tl ) ; 

3320nexti 
3330pri nt 
3340pri nt 

3350pri nt"Equi 1 ibri urn  mole  fractions  in  the  order  CH4,  CO,  C02,  H2,  H20  are:" 
3360fori=l to5 

3370  i fd$="£"ord$=" L"then3400 
3380  printusing".####  ",x(i); 

3390  goto341 0 
3400  printx(i); 

341 Onexti 
3420pri nt 
3430pri nt 

3440 i fd$=" 1 "ord$="L"then3470 

3450printusing"The  number  of  moles  of  gas  = #.####  per  mole  of  original  gas",g(t 

1 )/g 

3460goto3480 

3470pri nt"The  number  of  moles  of  gas  ="g(tl )/g;"per  mole  of  original  gas" 

3480 i fc$="p"orc$="P" then3510 

3490pri nt"The  pressure  ="g(tl )*p0/g;"atm  = ". 101325*g(tl )*p0/g; 

"MPa  per  mole  of  original  gas" 

3500goto3520 

3510print"The  volume  ="g(tl)*v0/g;"  liters  per  mole  of  original  gas" 

3520pri nt 
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3530oncgoto3540 , 3540,3590 , 3630 
3540  i fd$=ll£"ord$="  L"then3570 

3550pri ntusi ng"Sol id  carbon  is  present  in  the  amount  of  #.####  moles  per  mole  of 
original  gas",fnc/g 
3560goto3580 

3570print"Sol id  carbon  is  present  in  the  amount  of  "fnc/g;"moles  per  mole  of  original 

gas" 

3580 i fCl=2then3630 
3 5 9 0 i fd$="£"  or  d$="L"then3620 

3600printusing"Liquid  water  is  present  in  the  amount  of  #.####  moles  per  mole  of 
original  gas" ,n(8,tl )/g(0) 

361 0goto3630 

3620print"Liquid  water  is  present  in  the  amount  of  "n(8,tl )/g;"moles 
per  mole  of  original  gas" 

3630pri nt 

3640ifabs(k4-kl)/kl>le-4then3710 
3650ifabs(k5-k2)/k2>le-4then3710 
3660onc 1 go to3670 , 3670 , 3690 , 3690 
3670ifabs(k6-k3)/k3>le-4then3710 
3680goto3750 

3690if(k3-k3)/k3>le-4then3710 

3700goto3750 

37 1 Opr i nt"Equi 1 i bri um  constants  do  not  check  out" 

3720print"First  equilibrium  constant  ="kl;"vs."k4 

3730pri nt"Second  equilibrium  constant  ="k2;"vs."k5 

3740pri nt"Thi rd  equilibrium  constant  =" k3 ; "vs. 11  k6 

3 7 5 0 i ft>647. 3then3840 

3760onclgoto3770, 3790, 3770, 3790 

3 7 7 0 i fabs(dl )>5e-5then3800 

3780goto3850 

3790 i fdl <5e-5then3850 

3800print"Relationship  of  H20  partial  pressure  and  vapor  pressure  is  wrong" 

3810pri nt"Vapor  pressure  of  water  ="exp(£(7)) 

3820print"Partial  pressure  of  water  ="p*x(5) 

3830goto3850 

3840pri nt"Above  the  critical  temperature  for  water" 

3850pri nt 
3860' end 
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interactive  computer  program  written  in  the  BASIC  language  and  named  COLGAS.  The 
aim  of  this  work  was  to  provide  people  who  test  materials  in  coal-gasification-like 
atmospheres  an  easy  way  to  obtain  the  equilibrium  composition  of  their  gas  mixtures. 
A knowledge  of  computer  programming  is  not  required  in  order  to  use  the  program. 

A listing  of  the  program  is  given  and  also  six  sample  computer  calculations.  The 
phase  rule  is  applied  to  the  C-H-0  system  and  two  ternary  diagrams  are  shown 
illustrating  the  condensation  of  solid  carbon  and  liquid  water. 
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